Stabilization of ultra-short pulses in cubic nonlinear media 
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We study the propagation of ultra-short pulses in a cubic nonlinear medium. Using multiple-scale 
technique, we derive a new wave equation that preserves the nonlocal dispersion terms present in 
Maxwell's equations. As a result, we are able to understand how ultra-short nonlinear shocks are 
stabilized by dispersive terms. A delicate balance between dispersion and nonlinearity leads to a 
new type of solitary waves. Their stability is confirmed by numerical simulations of full Maxwell's 
equations. 

PACS numbers: 42.81.-i,42.25.Dd 



I. INTRODUCTION 

Recent progress in nonlinear optical systems has led to 
increasing development and use of ultra-short technolo- 
gies. Driven by the advent of novel detection techniques 
such as frequency resolved optical gating (FROG), the 
precise measurement of ultra-short pulses is now possi- 
ble. Understanding the dynamics of ultra-short pulses 
has been the subject of intensive research due to their po- 
tential uses for various fields, e.g., medical applications 
and imaging continuum recompression and control 
from highly nonlinear photonic crystal fibers |2|, L3 1 and 
the next generation of telecommunication systems. 

The cubic nonlinear Schrodinger equation (NLSE) has 
been the governing model for pulse propagation in non- 
linear medium. The NLSE, however, was shown to be in- 
adequate for describing ultra-short pulses. As the pulse 
power increases, i.e., the pulse becomes shorter, the de- 
lay in the response of the media to the excitation of 
the electric field starts to play a dominant role. The 
characteristics of the delay are not instantaneous, and 
hence they need to be described by convolution terms 
in governing equations. In contrast to differential opera- 
tors, which are local operators, convolution terms involve 
nonlocal operators and hence the resulting system is an 
integro-differential equation. The NLSE is a local equa- 
tion, namely, it only involves local terms approximating 
the delayed response. Thus, the NLSE remains valid only 
while the basic assumption, the existence of fairly broad 
pulses, is guaranteed Q|. 

Current optical technology allows to design state-of- 
the-art optical devices such as photonic crystals or mi- 
crostructured fibers whose structures are more complex 
than standard optical fibers. Experimental observations 
have shown that these new devices provide remarkable 
phenomena never seen in standard optical fibers. Some 
of the most prominent phenomena include creation of 
photonic band gaps, supercontinuum generation with 
relatively low pulse intensities, ultrashort pulse gener- 
atiompulse fission reminiscent of multisoliton breakup 
LJ- lil H , and the simultaneous third harmonic generation 
[g. While photonic crystal technology can successfully 



demonstrate the remarkable performance of the devices, 
very little is known in theoretical study. It is evident that 
one must develop a reliable mathematical model that can 
explain these new optical phenomena. In fact, any en- 
velope approximation that eventually leads to the cubic 
NLSE fails to exhibit the optical mechanism which evi- 
dently goes beyond the local soliton dynamics. 

A considerable amount of research efforts have been 
made to capture the complex optical phenomena that 
cannot be explained by the cubic NLSE. For short pulse 
dynamics, for example, one needs to introduce higher or- 
der dispersive terms, which lead to modified NLSE whose 
structure, although much more complex, still remains 
local. For certain materials, e.g. bulk silica, one can 
derive a local equation which describes the evolution of 
ultra-short pulses by making use of the specific form of 
the material susceptibility, 9|- Recently, the short 

pulse equation (SPE) was independently derived, which 
was proven to be a better approximation as the NLSE 
loses accuracy 10] . The SPE, however, was not shown 
to possess soliton solutions 0. Moreover, the SPE is not 
valid for all materials since the equation is derived tak- 
ing into account a particular type of the susceptibility of 
bulk silica. In this paper, considering general forms of 
susceptibility, we derive a new model without imposing 
any locality assumption and therefore, avoiding poten- 
tial loss of smoothness that could have helped to obtain 
stable solutions in SPE. While soliton theory has been 
developed in the context of locality, it is now our main 
task to investigate if nonlocal equations can provide sim- 
ilar phenomena. Indeed, we show that our new nonlocal 
equation also possesses solitary waves. Furthermore, we 
perform numerical simulations to confirm the stability of 
the solutions, which suggests that the soliton theory can 
be extended to systems with convolution terms. 

The paper is organized as follows. In section [HJ we 
present the derivation of the nonlinear nonlocal wave 
equation from Maxwell's equations. In section II I II us- 
ing multiple-scale technique, we derive a weakly nonlin- 
ear approximation to the wave equation introduced in 
section [H] In section IIVI we show that the absence of 
dispersion leads to the shock formation. In section we 
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further stabilize the optical shocks by introducing a dis- 
persive convolution term, and show our model possesses 
solitary waves. In section I VII we confirm the stability 
of the solitary wave solutions by numerical simulations. 
Finally, in section IVIII we summarize the main results of 
our analysis. 



II. DERIVATION OF NONLINEAR NONLOCAL 
WAVE EQUATION 

Prominent examples of effects that lead to convolution 
terms are optical phenomena due to the retardation in 
the response of the media to excitation by an oscillating 
electric field ^lj . Recently, the relation between nonlocal 
terms and quadratic solitons in quadratic nonlinear x^~ 
materials has been studied |l2(. Here, we consider cubic 
nonlinear materials and investigate their impacts on 
wave propagation. We start with Maxwell's equations: 



V x E 



dB 



VxH=j 



~&t' 



V • D = p, 

V-B = 0. (2.1) 



The fields E, B, D, H are real-valued vectors depending 
on the three space variables x, y, z and the time t; p is 
the charge density and j is the current density vector. We 
write $1, $2, $3 to denote the the first, second, and the 
third component of the vector The main results of this 
paper do not depend on the dimension of the considered 
system. Thus, for simplicity, the analysis is carried out 
in an one-dimensional context. This allows us to reduce 



J 



the Maxwell's equations to one dimensional system of 
equations by setting 



E = E 3 (x, t)e 3 , H = H 2 (x, t)e 2 . 



(2.2) 



In dielectric media, there exist neither free charges nor 
field sources, and hence we can set j = and p = 0. We 
also assume B = /ioH where po is a constant, which is 
appropriate for bulk silica. The relation between D and 
E is given by 



D, 



£qE 3 



(2.3) 



where eo is the permeability and P 3 will be defined later. 
These assumptions reduce the system of equations (|2.1f) 
to two one-dimensional equations: 



0H 2 dD 3 



dE 3 



dx dt ' dx ^° dt 
Using the above relations, we find the equation for E 3 , 



dH 2 



(2.4) 



d 2 E 3 8 2 E 3 d 2 P 3 



dx 2 



eo Mo- 



at 2 



Mo- 



di 2 



(2.5) 



Now we set c 2 = (/io £ o) 1 an d obtain the standard form 
of the wave equation: 



2 E 3 1 8 2 E 3 



dx 2 



dt 2 



c^di 2 



(2.6) 



Since the response of the medium is not instantaneous, 
we account for this retarded response in terms of convo- 
lution integrals and obtain 



Ps = e I X (1) {t-r)E 3 {x,T)dT + e / X (3) (* - n, t - t 2 , t - t 3 )E 3 (x, t x )E 3 {x, t 2 )E 3 (x, T 3 )d Tl dT 2 dT 3 . (2.7) 

I 



Note that this is where convolution terms enter naturally. (replacing E 3 by E and setting c = 1) 
We finally find one dimensional nonlinear wave equation 

I 



d 2 E d 2 E d 2 



dx 2 dt 2 dt 2 



X (1) {t-T)E{x,T)dT + J X {Z) {t-T 1 ,t~T 2 ,t~T 3 )E(x i T 1 )E{x i T 2 )E{x,T 3 )dT 1 dT 2 dT 3 ^j (2. 



Here, convolution terms are present in both the disper- 
sive and the nonlinear terms. Generally, the convolution 
integral in the dispersion has more significant impact on 
pulse propagation. The Fourier transform helps to con- 
vert the convolution present in the dispersive term to 
a product, however, even if the nonlinearity was local, 



it would lead to a convolution of the nonlinear term in 
Fourier domain. 
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III. NONLOCAL TERMS IN MULTI-SCALE 
EXPANSIONS 

Multiple-scale technique allows to derive asymptotic 
expansions for ordinary and partial differential equations 
|13j . The basic idea is to use the slow scales to remove 
resonances which are responsible for short validity range 
of the chosen expansion. In this paper, we focus on the 
case where dispersion and nonlinearity effects are weak 
and of the same order e. The method presented here is 
very general and hence, can be applied to other cases. In 



the presence of weak dispersion and nonlinearity, the cal- 
culations are more straightforward to perform since the 
zero-order operator is the unperturbed linear wave equa- 
tion whose kernel can be constructed by simple exponen- 
tial functions. This makes it easy to apply the Fredholm 
alternative leading to the solvability conditions which re- 
sult in the equations for the slow evolution of the system. 

Let us write the basic wave equation again, this time 
with a small dispersive and small nonlinear term on the 
right hand side: 



d 2 E 
dx 2 



8 2 E cP_ 



X {1) (t-T)E(x,T)dr 



X (3) (t ~Tl,t- T 2 ,t- T 3 )E{x,T 1 )E{x,T 2 )E(x,T 3 )dT 1 dT 2 dT 3 (3.1) 



r 



We introduce multiple scales only in the evolution vari- 
able z, and thus no additional scales in t. Then, in con- 
trast to the derivation of the NLSE, our analysis is valid 
even in the case where the Fourier transform is not con- 
centrated around one special frequency. This is of par- 
ticular importance for the study of very short pulses or 
phenomena involving frequency-mixing. The first step is 
to rewrite l|3.1fl in Fourier domain 

(J^ + <Ae = -eLU 2 X U(u,)E(x, u) ~ e-^M{E), 

(3.2) 

where we have introduced the nonlinear (and nonlocal) 
operator 

x,uj 1 )E(x,uj 2 )E(x,uj 3 ) 

x S(u> — loi — uj 2 — uj 3 ) du>i dus 2 dui 3 . (3-3) 
Now we perform a multiple scale expansion 

E{x,u) = E (x ,Xi,u)) + eEi(x ,xi,u)) H , 

xo = x, x\= ex, (3-4) 

and solve Eq. 113.211 order by order in e. The zero-order 
in e yields 



dx 2 



LU 2 ) Eq = 



(3.5) 



with the solution 

E {x ,x u u)=A {x 1 ,u)e i " x + B {x 1 ,u)e- i " x (3.6) 

where the term containing Ao corresponds to a wave mov- 
ing to the right and Bo to a wave moving to the left. 

Collecting the first order terms in e, we find the equa- 
tion for Ei as 
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d_ d_ 

dxo dxi 



— +u 2 )E l = -2— — E -lu 2 x (L) (cj)E 



(2k) 



jM(E ). 



(3.7) 



The question is now how to formulate the solvability con- 
dition for H3.7J) taking into account the nonlocal driving 
term on the right hand side of the equation. The fun- 
damental solutions of the linear operator on the l.h.s. of 
Eq. H3.7fl are again {exp(iwa;o), exp(— ituxo)}. Follow- 
ing the technique of multiple scales, we apply Fredholm 
alternative for removing possible resonances which arise 
from the terms in r.h.s. of Eq. I|3.7[) proportional to 
those fundamental solutions. Using the expression (|3.6|l . 
the triple product (omitting the arguments xq and x\ for 
simplicity) Eo(wi)Eq(u! 2 )Eo(lu — lj\ — oj 2 ) will yield eight 
terms. In each of those terms, we find all the possible 
combinations of u>\ and lj 2 that lead to the solution E\ 
proportional to either exp(itt>a;o) or exp(— iloxq). Those 
combinations are responsible for resonance and hence, 
must be removed. From Eq. (|3.(i|l . we find (the primes 
denote derivatives with respect to the slow variable x\) 



d d - 

■2q Q E (X0,XI,LU) 

oxq ox i 



-2i 



fa 



Thus, the contributions proportional to exp(za;a;o), 
exp(— iioxo) will appear in the equation for A' and B' , 
respectively. Let us now observe the first of the eight 
terms (omitting the evolution variable xq for the sake of 
easy typesetting) 

i (wi)i (cj 2 )i (^ -ojx- uj 2 )e luxa . 

The oscillations arise from e luJX ° and thus, this term will 
appear only in the equation for A' Q . The next term 

io(wi)io(w2)£o(w - wi - uj2)e^- u+2ui+2 ^ X0 

incorporates oscillations with different frequencies cor- 
responding to e l (-"+2^i+2w 2 )a;o Contributions to reso- 
nance terms occur if 

1. — lu + 2lui + 2lu 2 = lu responsible for contributions 
to the equation for A' . 
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2. — uj + 2lui + 2u>2 — —<-<J responsible for contributions 
to the equation for B' . 

Therefore, the solvability conditions lead to the system 
of equations for A' Q and B' Q , 

- 2iloA'^) = lo 2 x (1} A (lu) + ^M(A ) + 

35 (0) / x (3) (wi,w-w 2 ,0)i (wi)io(w-wi)^i 



(2tt)' 



(3.9) 



2iojB' (cj) = ^ 2 x (1) HA>M + ^^M(B )+ 



{2*f 



3i (0) / X (3) K^- LU2,0)B (uJl)Bo{uJ - LU^dU! 



+ 3B (uj) J x (3) (^i)-wi>w)i (wi)i (-wi) dwiV 

(3.10) 

The resulting system of equations for Ao and Bq describes 
the slow evolution of Eq. (|2.8(l . These provide a signif- 
icant simplification of the original equation (|2.8|l . which 
are the first order equations in x\ and, if we take the 
inverse Fourier transform, first order in t as well. Note 
that the resulting equations preserve the nonlocal struc- 
ture, which was our main purpose. This pair of equa- 
tions can be seen as an analog to the NLSE valid for 
broad pulses (without the assumption of a signal centered 
around a specific carrier frequency) or the SPE valid for 
short pulses (without the specific assumption about the 
profile of linear and nonlinear susceptibility). The only 
assumption necessary to obtain the above system was 
the smallness of the dispersive and the nonlinear terms. 
Our derivation also suggests how to proceed in the case 
where the dispersive term is not small. In this case, the 
solutions of the linear problem will involve a more compli- 
cated dispersion relation which will give rise to modified 
resonance conditions for the nonlinearity. 

Note that Eqs. (|3.9|l . (|3.10() also account for coupling 
between the forward and backward moving waves Aq, 
Bq, respectively. The coupling is weak in the sense that 
it only consists of integrated quantities. In some cases, 
however, e.g. when the susceptibilities have an additional 
dependence on the evolution variable x such as in certain 
photonic crystal structures, the coupling between Aq and 
Bq can become important. 



IV. THE DISPERSION-FREE CASE 

In the search of soliton solutions of the underlying non- 
local equations, we first consider a particular case, ab- 
sence of linear dispersion x^ 1 ' = and constant nonlinear 
susceptibility x^ = 1- Then, we take the inverse Fourier 



transform of (|3.9() , i|3.10[l to obtain differential equations 
in time domain. 

We first notice that B' o (0, x\) — implies Bq(0) = 
Bq(cj — 0,xi) is constant as a function of x\ and the 
corresponding statement also holds for A (ui = 0, Xi). 
For convenience of notation, let us define 

2irB zcto = B (0), 2ttA zcio = i (0). 

After straightforward calculation we find that 

1 



A int := 



-Bint. := 



(2tt)2 

1 



Ao(u>i)A \-u>i) dcji, 
Bo(wi)B (-wi) dwi (4.1) 



are constants of motion as well. Therefore, after Fourier 
transform back in time domain we have (Aq = AqIxA) 
and Bq — Bq[x, t)) 

■t^-Aq + l^- (A 3 + 3B zcro A 2 + 3B int A Q ) = 0, (4.2) 
oxi 2 at 

■T- B ° -lirA B o+ 3A ^oB 2 + 3A int B ) = 0. (4.3) 
ox\ 2 ot 

As mentioned before, the forward and backward propa- 
gating waves Aq, Bq are coupled. This coupling, how- 
ever, is weak in a sense that it occurs only through 
the constants of motions. By choosing initial condition 
-Bo^i = 0, t) =0 we can eliminate the backward prop- 
agating wave Bq entirely at the leading order and the 
equation for A is reduced to a Burger's type equation 



9 A 1 9 A3 n 



with the (implicit) solution 



A (t, Xl ) = f [t--A 2 Q Xl 



(4.4) 



(4.5) 



Here, the function / is determined by the initial condi- 
tion of the pulse at t = 0. For typical pulse profiles (e.g., 
a Gaussian pulse), we see that shocks will be generated 
eventually. This was expected as we excluded all the lin- 
ear dispersion terms from our considerations that could 
have helped to smoothen the solution. 



V. SOLITARY WAVES 

In order to regularize the shocks and further obtain 
stable wave propagation, we bring back a small dispersion 
term (x^ is still set to be constant for simplicity) and 
find 



dx 2 



dt 2 



E(x,t) 



J X {1) (t-r)E(x,T)dT 



(5.1) 
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Now we apply the same analysis used before and obtain 
an additional nonlocal term in Eqs. (|4.2I) . 14. 3[) . 

^T^o + 2 ^ (^o + 3B zcj:o Aq + 3B int A 

+ J x {1 \t-T)A { Xl ,T)dT^ = 0, (5.2) 

t; — B — -— [B Q + 3^4 zcro -B + 3A ln tB 
oari 2 at 

+ J X {1) {t-r)B Q { Xl ,T)dT^ =0. (5.3) 

Again we assume the back-scattered part Bq is sup- 
pressed. Motivated by the recent developments e.g., in 
photonic crystal structures a wide range of susceptibili- 
ties can be engineered, we search for the profiles of Ao 
and x which can stabilize the shock formation. One 
simple approach to finding a stationary solution is to set 
dAo/dxi = 0. This implies we need to find Ao and x^ 1 -* 
such that 

Al + J x (1) {t~T)A ( Xl ,T)dT = 0. (5.4) 

Denoting x = F(x )> we rewrite Eq. 1)5. 4f> in Fourier 
domain, 



^)+ X (1) io = 0. 
If we assume a Lorentz profile for Aq 

A (t) 



i + (/3ty 



(5.5) 



(5.6) 



we can solve Eq. 1)5.4( 1 using Fourier transform. The 
Fourier transform of Aq is 



ioH = vr^e-H//? 
P 



(5.7) 



and, using the residue theorem, we obtain the Fourier 
transform of Aq, 



Most of the physical examples for Maxwell's equations, 
however, require a modulated initial condition. There- 
fore, we consider a pulse of the form 



A (t) = A (t) cos(w t) 



(5.10) 



which will have its maximal frequency coefficients in the 
neighborhood of ujq. We extend our analysis to this case 
of a modulated initial condition and find a susceptibility 
for which the given pulse will maintain its stable propa- 
gation. Since the Fourier transforms of Aq and Aq can 
be found directly as 

T{Aq) = - (A (v - cj ) + A (w + uo)^J (5.11) 

HAD = |(^)(w + 3wo)+aF(^)(w + wo) 

+ 3T(A 3 )(u -UQ+ T{J§)(w - 3o> )) (5.12) 

we find the susceptibility Hx^) for the modulated soli- 
tary wave as 



Hx {1) ) = 



HAD 
HM 



(5.13) 



It is important to note that the presented susceptibili- 
ties (I5.9|) and (|5.13|l do not correspond to susceptibili- 
ties of any known material. Although they need to be 
engineered artificially, recent progress in manufacturing 
particular photonic crystal structures 0, 0, 0] makes 
it reasonable to expect the possibility to design such and 
similar susceptibilities in a not too distant future. On 
the other hand, the above analysis offers a variety of dif- 
ferent approaches for regularizing shocks to obtain stable 
propagation of nonlinear waves. One possibility, for ex- 
ample, would be to introduce a dependence of the slow 
variable x\ to x^ ( as i n certain photonic crystal struc- 
tures). This increases the mathematical difficulty to find 
solutions to the system of equations 1)5. 2|) , 15.311 but leads 
to more freedom in the choice of Aq and x^ ■ 



HAD(") = ^e-M//V + 3/?M + 3/3 2 ). (5.8) 

We find that Eq. (|5.6() is a stationary solution provided 
that the susceptibility takes the following parabolic form: 



i^ 

"8/3 2 



{u? + 3/3M +3/3 2 



(5.9) 



We expect that a pulse of this particular shape, for the 
given susceptibility will not undergo changes at the lead- 
ing order, hence propagates as a stable nonlinear wave. 
The dispersion compensates for the effect of nonlocal 
nonlinear term not only in a small frequency range but 
for all possible values of cu. Due to this balance between 
convolution terms in the nonlinearity and dispersion, we 
can obtain a stable solitary wave. 



VI. NUMERICAL RESULTS 

In order to show the stability of the solitary wave 1)5.6)) 
numerically, we use its analytic expression together with 
the corresponding susceptibility l|5.9|l for simulations of 
the nonlocal nonlinear wave equation 12.8(1 . Here, we 
chose a — 0.2 and (3 = 0.75 leading to an effective e of 
around 0.07. Transmitting the solitary wave to a dis- 
tance x = 25, which is longer than C(l/e), enables us to 
capture all three important effects: Pulse distortion due 
to dispersion, pulse distortion due to nonlinearity, and 
stable propagation of the nonlinear soliton due to a non- 
local balance of dispersive and nonlinear effects . Figure 
^presents the pure linear and nonlinear effects. In Fig. |2 
we present the propagation of the solitary wave combin- 
ing the linear and nonlinear effects. Although the char- 
acteristics of the soliton are resulted from a first-order 



15 20 25 30 35 
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FIG. 1: Figure on top: The diamonds show the solution of 
Maxwell's wave equation 12,81 in the linear case x' 3 ' = 
versus the initial pulse (dashes) in a moving frame. The linear 
susceptibility x leads to dispersion. Bottom figure: The 
diamonds represents the solution of Maxwell's wave equation 
12.811 in the purely nonlinear case x^ = versus the initial 
pulse (dashes) in a moving frame. It clearly illustrates the 
formation of a shock. 



approximation, the numerical simulations show that the 
balance of nonlinear convolution terms and dispersion 
generates a stable solution whose shape changes only 
slightly during the propagation. 

The second numerical simulation shows the stable 
propagation of the modulated initial pulse given by <|5.10|) 
and corresponding susceptibility l|5.13l) . Here, we chose 
simulation parameters a — 0.2 and (3 = 1.0 and ujq = 5.0. 
The allowed frequency range was restricted to frequen- 
cies below 2 loq to suppress generation of third harmon- 
ics. Figure |31 shows how well nonlinearity and dispersion 
are balanced in this case. The pulse experiences only 
very small deviations from its original shape that are due 
to higher-order effects. From this figure it is also clear 
that this solitary wave represents an ultrashort pulse, 
hence does not fall into the regime of the cubic nonlinear 
Schrodinger equation. It is a new type of ultra short soli- 
tary wave that can exist due to the balance of dispersion 
and nonlinearity in a very broad frequency range. 



15 20 25 30 35 

t 

FIG. 2: Propagation of the solitary wave. The diamonds 
represent the solution of Maxwell's wave equation with both 
dispersive and nonlinear terms. The stability of the pulse 
propagation is due to the balance of both effects. 




20 22 24 26 28 30 

t 

FIG. 3: Propagation of the modulated solitary wave. The 
diamonds represent the solution of Maxwell's wave equation 
with both dispersive and nonlinear terms, the continuous line 
the initial condition in a moving frame. The dashed line is 
the linear solution of Maxwell's equations with the same ini- 
tial condition. The figure shows how nonlinearity allows the 
solitary wave to maintain its initial shape. 



VII. CONCLUSION 

We derived a slowly varying envelope equation for 
ultra-short pulses from the cubic nonlinear Maxwell's 
equation. The derivation was carried out without impos- 
ing any locality assumption unlike NLSE, which allows 
to preserve the phenomena that are due to the retarded 
response of the media in ultra-short pulse dynamics. We 
also did not assume any specific types of material suscep- 
tibility and hence, our new model can be applied more 
generally than the SPE. We showed that for certain sus- 
ceptibility of material the resulting approximative equa- 
tions possess a solitary wave solution due to a (nonlocal) 
balance between dispersion and nonlinearity. The stable 
propagation of this solitary wave is confirmed by numer- 
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